home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / usenet / sources / volume90 / aplictns / listplot / part03 < prev   
Internet Message Format  |  1990-08-23  |  29KB
  1. Path: abcfd20.larc.nasa.gov!amiga-request
  2. From: amiga-request@abcfd20.larc.nasa.gov (Amiga Sources/Binaries Moderator)
  3. Subject: v90i244: ListPlot - 2d plotting program, Part03/03
  4. Reply-To: Anthony M. Richardson <amr@dukee.egr.duke.edu>
  5. Newsgroups: comp.sources.amiga
  6. Message-ID: <comp.sources.amiga:v90i244@abcfd20.larc.nasa.gov>
  7. Date: 23 Aug 90 01:03:49 GMT
  8. Approved: tadguy@uunet.UU.NET (Tad Guy)
  9. X-Mail-Submissions-To: amiga@uunet.uu.net
  10. X-Post-Discussions-To: comp.sys.amiga
  11.  
  12. Submitted-by: Anthony M. Richardson <amr@dukee.egr.duke.edu>
  13. Posting-number: Volume 90, Issue 244
  14. Archive-name: applications/listplot/part03
  15.  
  16. #!/bin/sh
  17. # This is a shell archive.  Remove anything before this line, then unpack
  18. # it by saving it into a file and typing "sh file".  To overwrite existing
  19. # files, type "sh file -c".  You can also feed this as standard input via
  20. # unshar, or by typing "sh <file", e.g..  If this archive is complete, you
  21. # will see the following message at the end:
  22. #        "End of archive 3 (of 3)."
  23. # Contents:  Csrc/plotting.c
  24. # Wrapped by tadguy@abcfd20 on Wed Aug 22 21:03:22 1990
  25. PATH=/bin:/usr/bin:/usr/ucb ; export PATH
  26. if test -f 'Csrc/plotting.c' -a "${1}" != "-c" ; then 
  27.   echo shar: Will not clobber existing file \"'Csrc/plotting.c'\"
  28. else
  29. echo shar: Extracting \"'Csrc/plotting.c'\" \(25452 characters\)
  30. sed "s/^X//" >'Csrc/plotting.c' <<'END_OF_FILE'
  31. X/*
  32. X * plotting.c
  33. X *
  34. X * plotting routines.
  35. X *
  36. X */
  37. X#include <stdio.h>
  38. X#include <errno.h>
  39. Xextern int errno;
  40. X#include <assert.h>
  41. X#include <string.h>
  42. X#include <ctype.h>
  43. X#include <signal.h>
  44. X#include "datatypes.h"
  45. X
  46. Xvoid    ErrorExit();
  47. X
  48. X#ifndef    min
  49. X#    define    min(A,B)    A<B? A : B
  50. X#endif
  51. X#ifndef    max
  52. X#    define    max(A,B)    A>B? A : B
  53. X#endif
  54. X
  55. Xextern bool    GraphicsInProgress;
  56. X
  57. X/* GoldenRatio = (1+Sqrt(5))/2    */
  58. X#define    GOLDENRATIO    1.6180339887499
  59. X#define ASPECTRATIO    1.0/GOLDENRATIO
  60. Xextern VALUE    V_AngularUnit;
  61. Xextern VALUE    V_AnnotationScale;
  62. Xextern VALUE    V_AspectRatio;
  63. Xextern VALUE    V_Boxed;
  64. Xextern VALUE    V_Domain;
  65. Xextern VALUE    V_Gridding;
  66. Xextern VALUE    V_LabelScale;
  67. Xextern VALUE    V_LineColor;
  68. Xextern VALUE    V_LineStyle;
  69. Xextern VALUE    V_Orientation;
  70. Xextern VALUE    V_Origin;
  71. Xextern VALUE    V_PlotColor;
  72. Xextern VALUE    V_PlotDevice;
  73. Xextern VALUE    V_PlotJoined;
  74. Xextern VALUE    V_PlotPoints;
  75. Xextern VALUE    V_PlotTitle;
  76. Xextern VALUE    V_PlotType;
  77. Xextern VALUE    V_PointScale;
  78. Xextern VALUE    V_PointSymbol;
  79. Xextern VALUE    V_PolarVariable;
  80. Xextern VALUE    V_Range;
  81. Xextern VALUE    V_SubPages;
  82. Xextern VALUE    V_SupplyAbscissa;
  83. Xextern VALUE    V_TitleScale;
  84. Xextern VALUE    V_UseInputFile;
  85. Xextern VALUE    V_Verbose;
  86. Xextern VALUE    V_ViewPort;
  87. Xextern VALUE    V_XLabel;
  88. Xextern VALUE    V_XTick;
  89. Xextern VALUE    V_YLabel;
  90. Xextern VALUE    V_YTick;
  91. X
  92. Xchar    *DeviceTypes[] = {
  93. X    "amiga",
  94. X    "printer",
  95. X    "iff",
  96. X    "hp",
  97. X    "aegis",
  98. X    "postscript",
  99. X    (char *)NULL
  100. X};
  101. X
  102. Xtypedef struct line_style {
  103. X        int    ls_type;
  104. X        int    ls_nelms;
  105. X        int    *ls_space;
  106. X        int    *ls_mark;
  107. X    } LINESTYLE;
  108. XLIST    LineStyleList;
  109. XLINESTYLE    *LineStyles; int NLineStyles;
  110. Xint        *PointCodes, NPointCodes;
  111. Xint    space0 = 0, mark0 = 0, space1 = 1500, mark1 = 1500;
  112. X
  113. XRECT    UserRect;
  114. XRECT    ViewPort = { 0.1,0.1, 0.9,0.9 };
  115. XRECT    WorldRect = {0.0,0.0, 1.0,1.0};
  116. X
  117. X
  118. Xvoid
  119. XListPlot(Data, NTuples, TupleSize)
  120. XFLOAT    **Data;
  121. Xint    NTuples, TupleSize;
  122. X{
  123. X    void InitializeDevice();
  124. X    void InitializePlottingEnv();
  125. X#ifdef    ANSI_C
  126. X    void Plot2D(FLOAT **, int , int );
  127. X    void PolarPlot(FLOAT **, int , int );
  128. X#else
  129. X    void Plot2D();
  130. X    void PolarPlot();
  131. X#endif
  132. X
  133. X    InitializeDevice();
  134. X    GraphicsInProgress = TRUE;
  135. X
  136. X    InitializePlottingEnv(Data, NTuples, TupleSize);
  137. X
  138. X    if (Vstrcmp(V_PlotType,"polar")) {
  139. X        PolarPlot(Data, NTuples, TupleSize);
  140. X    } else {
  141. X        Plot2D(Data, NTuples, TupleSize);
  142. X    }
  143. X    plend();
  144. X}
  145. X
  146. X
  147. Xvoid
  148. XInitializeDevice()
  149. X{
  150. X    register int    i;
  151. X    register char    *DeviceStr;
  152. X
  153. X
  154. X    plorient(Vstrcmp(V_Orientation, "portrait")? 1 : 0);
  155. X    plselect(stdout);    /* write to stdout */
  156. X        
  157. X
  158. X    DeviceStr = V_PlotDevice.val_u.udt_string;
  159. X    for (i=0; DeviceTypes[i]; i++) 
  160. X        if (strcmp(DeviceStr, DeviceTypes[i]) == 0) {
  161. X            plbeg(
  162. X                i+1,    /* device code    */
  163. X                (int)V_SubPages.val_u.udt_interval.int_lo, /* nx */
  164. X                (int)V_SubPages.val_u.udt_interval.int_hi /* ny */
  165. X            );
  166. X            return;
  167. X        }
  168. X    fprintf(stderr,"(InitializeDevice) Unknown device type: %s\n", DeviceStr);
  169. X    ErrorExit();
  170. X}
  171. X
  172. Xvoid
  173. XInitializePlottingEnv(Data, NTuples, TupleSize)
  174. XFLOAT    **Data;
  175. Xint    NTuples;
  176. Xint    TupleSize;
  177. X{
  178. X    int    i, j;
  179. X    FLOAT    *Ptr;
  180. X    void    InitLineStyles();
  181. X    void    InitPointCodes();
  182. X
  183. X    /* viewport    */
  184. X    ViewPort = VtoRect(V_ViewPort);
  185. X
  186. X    /* user window        */
  187. X    UserRect.XMin = MAX_FLOAT;
  188. X    UserRect.YMin = MAX_FLOAT;
  189. X    UserRect.XMax = -MAX_FLOAT;
  190. X    UserRect.YMax = -MAX_FLOAT;
  191. X    
  192. X    if (VtoBoolean(V_SupplyAbscissa)) {
  193. X        UserRect.XMin = 0.0;
  194. X        UserRect.XMax = (FLOAT)(NTuples-1);
  195. X    } else {
  196. X        for (i=0; i<NTuples; i++) {
  197. X            if (Data[i][0] < UserRect.XMin)
  198. X                UserRect.XMin = Data[i][0];
  199. X            if (Data[i][0] > UserRect.XMax)
  200. X                UserRect.XMax = Data[i][0];
  201. X        }
  202. X    }
  203. X    for (i=0; i<NTuples; i++) 
  204. X        for (j=TupleSize-1,Ptr= &Data[i][1]; j>0; --j,Ptr++) {
  205. X            if (*Ptr < UserRect.YMin)
  206. X                UserRect.YMin = *Ptr;
  207. X            if (*Ptr > UserRect.YMax)
  208. X                UserRect.YMax = *Ptr;
  209. X        }
  210. X        
  211. X    
  212. X    /* font scaling        */
  213. X
  214. X    /* plot type        */
  215. X
  216. X    /* grid            */
  217. X
  218. X    /* line styles        */
  219. X    InitLineStyles();
  220. X    
  221. X
  222. X    /* line colors        */
  223. X
  224. X    /* point symbols    */
  225. X    InitPointCodes();
  226. X
  227. X    /* window type        */
  228. X
  229. X    /* X label        */
  230. X
  231. X    /* Y label        */
  232. X
  233. X    /* Plot title        */
  234. X
  235. X    /* aspect ratio        */
  236. X}
  237. X
  238. X#define    isMark(C)    (((C) == 'M' || (C) == 'm')? TRUE : FALSE)
  239. X#define    isSpace(C)    (((C) == 'S' || (C) == 's')? TRUE : FALSE)
  240. X#define    MarkLength(C)    (C) == 'M'? 1000 : 250
  241. X#define    SpaceLength(C)    (C) == 'S'? 1000 : 250
  242. X#define    MAX_MARKS    32
  243. X
  244. Xvoid
  245. XDecodeLineStyle(StyleStr, Mark, Space, N)
  246. Xchar    *StyleStr;
  247. Xint    **Mark;
  248. Xint    **Space;
  249. Xint    *N;
  250. X{
  251. X    register char    *S;
  252. X    int    macc, sacc;
  253. X    int    i, j, k;
  254. X    int    mark[MAX_MARKS];
  255. X    int    space[MAX_MARKS];
  256. X    bool    ProcessingMark;
  257. X
  258. X    /*
  259. X     * 'S' = 1000   micron space
  260. X     * 's' =  250   micron space
  261. X     * 'M' = 1000     micron line
  262. X     * 'm' =  250    micron line
  263. X     */
  264. X    ProcessingMark = TRUE; /* to satisfy some compilers */
  265. X    if (isMark(*StyleStr))
  266. X        ProcessingMark = TRUE;
  267. X    else if (isSpace(*StyleStr))
  268. X        ProcessingMark = FALSE;
  269. X    else {
  270. X        fprintf(stderr,"(DecodeLineStyle) Invalid Line style: \"%s\"\n", StyleStr);
  271. X        ErrorExit();
  272. X    }
  273. X    memset((char *)mark, 0, MAX_MARKS * sizeof(int));
  274. X    memset((char *)space, 0, MAX_MARKS * sizeof(int));
  275. X    for (macc=0,sacc=0,i=0,j=0,S=StyleStr ; *S && i<MAX_MARKS && j<MAX_MARKS; ++S) {
  276. X        if (isMark(*S) && ProcessingMark) {
  277. X            macc += MarkLength(*S);
  278. X        } else if (isMark(*S) && !ProcessingMark) {
  279. X            space[j++] = sacc;
  280. X            macc = MarkLength(*S);
  281. X            ProcessingMark = TRUE;
  282. X        } else if (isSpace(*S) && ProcessingMark) {
  283. X            mark[i++] = macc;
  284. X            sacc = SpaceLength(*S);
  285. X            ProcessingMark = FALSE;
  286. X        } else if (isSpace(*S) && !ProcessingMark) {
  287. X            sacc += SpaceLength(*S);
  288. X        } else if (isspace(*S)) {
  289. X            continue;
  290. X        } else {
  291. X            fprintf(stderr,"(DecodeLineStyle) Invalid line style code '%c' in \"%s\"\n", *S, StyleStr);
  292. X            ErrorExit();
  293. X        }
  294. X    }
  295. X    if (ProcessingMark)
  296. X        mark[i++] = macc;
  297. X    else
  298. X        space[j++] = sacc;
  299. X
  300. X    k = max(i,j);
  301. X
  302. X    if (((*Mark) = (int *)calloc(k, sizeof(int))) == (int *)NULL) {
  303. X        perror("(DecodeLineStyle) ");
  304. X        ErrorExit();
  305. X    }
  306. X    if (((*Space) = (int *)calloc(k, sizeof(int))) == (int *)NULL) {
  307. X        perror("(DecodeLineStyle) ");
  308. X        ErrorExit();
  309. X    }
  310. X    for (i=0; i<k; i++) {
  311. X        (*Mark)[i] = mark[i];
  312. X        (*Space)[i] = space[i];
  313. X    }
  314. X    *N = k;
  315. X}
  316. X
  317. Xvoid
  318. XInitLineStyles()
  319. X{
  320. X    int        i;
  321. X    LATOM        *A;
  322. X    LINESTYLE    *LS;
  323. X    int        *Mark, *Space;
  324. X    int        N;
  325. X#ifdef    ANSI_C
  326. X    LINESTYLE    *LSAlloc(int , int *, int *);
  327. X#else
  328. X    LINESTYLE    *LSAlloc();
  329. X#endif
  330. X
  331. X    LineStyleList.list_n = 0;
  332. X    LineStyleList.list_head = (LATOM *)NULL;
  333. X    LineStyleList.list_tail = (LATOM *)NULL;
  334. X
  335. X    LS = LSAlloc(0, (int *)NULL, (int *)NULL);
  336. X    Append(&LineStyleList, (generic *)LS);
  337. X
  338. X    for (i=1,A=V_LineStyle.val_u.udt_set.list_head; A; A=A->la_next,i++) {
  339. X        DecodeLineStyle(
  340. X            (char *)A->la_p,
  341. X            &Mark,
  342. X            &Space,
  343. X            &N
  344. X        );
  345. X        LS = LSAlloc(N, Mark, Space);
  346. X        Append(&LineStyleList, (generic *)LS);
  347. X    }
  348. X
  349. X    NLineStyles = i;
  350. X    if ((LineStyles = (LINESTYLE *)calloc(NLineStyles, sizeof(LINESTYLE))) == (LINESTYLE *)NULL) {
  351. X        perror("(InitLineStyles) ");
  352. X        ErrorExit();
  353. X    }
  354. X    for (i=0,A=LineStyleList.list_head; A; A=A->la_next,i++) {
  355. X        LineStyles[i].ls_mark = ((LINESTYLE *)(A->la_p))->ls_mark;
  356. X        LineStyles[i].ls_space = ((LINESTYLE *)(A->la_p))->ls_space;
  357. X        LineStyles[i].ls_nelms = ((LINESTYLE *)(A->la_p))->ls_nelms;
  358. X    }
  359. X}
  360. X
  361. X
  362. XLINESTYLE    *
  363. XLSAlloc(n, Mark, Space)
  364. Xint    n;
  365. Xint    *Mark, *Space;
  366. X{
  367. X    register LINESTYLE    *LS;
  368. X
  369. X    if ((LS = (LINESTYLE *)calloc(1, sizeof(LINESTYLE))) == (LINESTYLE *)NULL) {
  370. X        perror("(LSAlloc) ");
  371. X        assert(LS != (LINESTYLE *)NULL);
  372. X        exit(errno);
  373. X    }
  374. X    LS->ls_mark = Mark;
  375. X    LS->ls_space = Space;
  376. X    LS->ls_nelms = n;
  377. X    return(LS);
  378. X}
  379. X
  380. Xvoid
  381. XInitPointCodes()
  382. X{
  383. X    register LATOM    *A;
  384. X    register int    i;
  385. X    extern int    *PointCodes, NPointCodes;
  386. X    extern VALUE    V_PointSymbol;
  387. X
  388. X    if (!isSet(V_PointSymbol)) {
  389. X        NPointCodes = 0;
  390. X        return;
  391. X    }
  392. X
  393. X    /* count # codes    */
  394. X    for (i=0,A=V_PointSymbol.val_u.udt_set.list_head; A; A=A->la_next, i++)
  395. X        ;
  396. X
  397. X    NPointCodes = i;
  398. X    if (NPointCodes == 0)
  399. X        return;
  400. X
  401. X    if ((PointCodes = (int *)calloc(NPointCodes, sizeof(int))) == (int *)NULL) {
  402. X        perror("(InitPointCodes) ");
  403. X        assert(PointCodes != (int *)NULL);
  404. X        ErrorExit();
  405. X    }
  406. X    for (i=0,A=V_PointSymbol.val_u.udt_set.list_head; A; A=A->la_next, i++) {
  407. X        PointCodes[i] = atoi((char *)A->la_p);
  408. X    }
  409. X}
  410. X
  411. X
  412. Xvoid
  413. XPlot2D(Data, NTuples, TupleSize)
  414. XFLOAT    **Data;
  415. Xint    NTuples;
  416. Xint    TupleSize;
  417. X{
  418. X    int    i, j;
  419. X    FLOAT    *X, *Y;
  420. X    FLOAT    x, y;
  421. X    FLOAT    x0, y0;
  422. X    FLOAT    xtick, ytick;
  423. X    FLOAT    Xmin,Xmax, Ymin,Ymax;
  424. X    FLOAT    MedianX, StdDevX;
  425. X    FLOAT    MedianY, StdDevY;
  426. X    int    nxsub, nysub;
  427. X    bool    AutoRange, AutoDomain;
  428. X    bool    LogX, LogY;
  429. X    char    Xopt[16], Yopt[16];
  430. X#ifdef    ANSI_C
  431. X    void    StatisticsOfX(
  432. X            FLOAT **,
  433. X            int,
  434. X            int ,
  435. X            FLOAT *,
  436. X            FLOAT *
  437. X        );
  438. X    void    StatisticsOfY(
  439. X            FLOAT **,
  440. X            int,
  441. X            int ,
  442. X            FLOAT *,
  443. X            FLOAT *
  444. X        );
  445. X    double    log10(double);
  446. X#else
  447. X    void    StatisticsOfX();
  448. X    void    StatisticsOfY();
  449. X    double    log10();
  450. X#endif
  451. X
  452. X    pladv(0);
  453. X    ViewPort = VtoRect(V_ViewPort);
  454. X
  455. X    if (isString(V_AspectRatio)) {
  456. X        /* automatic --> STRETCH */
  457. X        plvpor(ViewPort.XMin, ViewPort.XMax, ViewPort.YMin, ViewPort.YMax);
  458. X/*debug
  459. X        plvsta();
  460. Xdebug*/
  461. X    } else {
  462. X        assert(isDbl(V_AspectRatio));
  463. X        plvasp(VtoDbl(V_AspectRatio));
  464. X    }
  465. X    
  466. X    LogX = (Vstrcmp(V_PlotType, "loglin") || Vstrcmp(V_PlotType,"loglog"))?
  467. X            TRUE : FALSE;
  468. X    LogY = (Vstrcmp(V_PlotType, "linlog") || Vstrcmp(V_PlotType,"loglog"))?
  469. X            TRUE : FALSE;
  470. X
  471. X    /* window bounds */
  472. X    AutoDomain = FALSE;
  473. X    if (isString(V_Domain) && Vstrcmp(V_Domain, "All")) {
  474. X        /* all points    */
  475. X        Xmin = UserRect.XMin;
  476. X        Xmax = UserRect.XMax;
  477. X    } else if (isString(V_Domain) && Vstrcmp(V_Domain,"Automatic")) {
  478. X        /* automatic     */
  479. X        AutoDomain = TRUE;
  480. X        StatisticsOfX(Data, NTuples, TupleSize, &MedianX, &StdDevX);
  481. X        x = MedianX - 3.0*StdDevX;
  482. X        Xmin = UserRect.XMin < x? x : UserRect.XMin;
  483. X        x = MedianX + 3.0*StdDevX;
  484. X        Xmax = UserRect.XMax > x? x : UserRect.XMax;
  485. X    } else {
  486. X        assert(isInterval(V_Domain));
  487. X        Xmin = V_Domain.val_u.udt_interval.int_lo;
  488. X        Xmax = V_Domain.val_u.udt_interval.int_hi;
  489. X    }
  490. X    if (LogX) {
  491. X        if (Xmin == 0.0 || Xmax == 0.0) {
  492. X            /* error */
  493. X            fprintf(stderr,"ListPlot: Trying to find log10(0) in data! Stop.\n");
  494. X            ErrorExit();
  495. X        }
  496. X        Xmin = log10((double)Xmin);
  497. X        Xmax = log10((double)Xmax);
  498. X    }
  499. X
  500. X    AutoRange = FALSE;
  501. X    if (isString(V_Range) && Vstrcmp(V_Range, "All")) {
  502. X        /* all */
  503. X        Ymin = UserRect.YMin;
  504. X        Ymax = UserRect.YMax;
  505. X    } else if (isString(V_Range) && Vstrcmp(V_Range,"Automatic")) {
  506. X        /* automatic     */
  507. X        AutoRange = TRUE;
  508. X        StatisticsOfY(Data, NTuples, TupleSize, &MedianY, &StdDevY);
  509. X        y = MedianY - 3.0*StdDevY;
  510. X        Ymin = UserRect.YMin < y? y : UserRect.YMin;
  511. X        y = MedianY + 3.0*StdDevY;
  512. X        Ymax = UserRect.YMax > y? y : UserRect.YMax;
  513. X    } else {
  514. X        assert(isInterval(V_Range));
  515. X        Ymin = V_Range.val_u.udt_interval.int_lo;
  516. X        Ymax = V_Range.val_u.udt_interval.int_hi;
  517. X    }
  518. X    if (LogY) {
  519. X        if (Ymin == 0.0 || Ymax == 0.0) {
  520. X            /* error */
  521. X            fprintf(stderr,"ListPlot: Trying to find log10(0) in data! Stop.\n");
  522. X            ErrorExit();
  523. X        }
  524. X        Ymin = log10((double)Ymin);
  525. X        Ymax = log10((double)Ymax);
  526. X    }
  527. X    plwind(
  528. X        Xmin,
  529. X        Xmax,
  530. X        Ymin,
  531. X        Ymax
  532. X    );
  533. X    plschr(0., VtoDbl(V_AnnotationScale));
  534. X
  535. X
  536. X    /* X and Y ticks */
  537. X    if (isString(V_XTick)) {
  538. X        /* automatic */
  539. X        xtick = 0.;
  540. X        nxsub = 0;
  541. X    } else {
  542. X        xtick = V_XTick.val_u.udt_interval.int_lo;
  543. X        if (LogX) xtick = log10((double)xtick);
  544. X        nxsub = V_XTick.val_u.udt_interval.int_hi;
  545. X    }
  546. X    if (isString(V_YTick)) {
  547. X        /* automatic */
  548. X        ytick = 0.;
  549. X        nysub = 0;
  550. X    } else {
  551. X        ytick = V_YTick.val_u.udt_interval.int_lo;
  552. X        if (LogY) ytick = log10((double)ytick);
  553. X        nysub = V_YTick.val_u.udt_interval.int_hi;
  554. X    }
  555. X
  556. X    /* Axis options    */
  557. X    Xopt[0] = (char)NULL;
  558. X    Yopt[0] = (char)NULL;
  559. X    if (isBoolean(V_Boxed) && VtoBoolean(V_Boxed)) {
  560. X        strcat(Xopt, "bc");
  561. X        strcat(Yopt, "bc");
  562. X    } else if (isString(V_Boxed)) {
  563. X        if (Vstrchr(V_Boxed, 'b'))
  564. X            strcat(Xopt, "b");
  565. X        if (Vstrchr(V_Boxed, 't'))
  566. X            strcat(Xopt, "c");
  567. X        if (Vstrchr(V_Boxed, 'r'))
  568. X            strcat(Yopt, "c");
  569. X        if (Vstrchr(V_Boxed, 'l'))
  570. X            strcat(Yopt, "b");
  571. X    } else {
  572. X        ;
  573. X    }
  574. X    if (LogX)
  575. X        strcat(Xopt, "l");
  576. X    if (LogY)
  577. X        strcat(Yopt, "l");
  578. X    strcat(Xopt, "nst");
  579. X    strcat(Yopt, "nstv");
  580. X
  581. X    if (isString(V_Origin) && Vstrcmp(V_Origin, "Median")) {
  582. X        if (AutoDomain == FALSE) {
  583. X            StatisticsOfX(Data,NTuples,TupleSize,&MedianX,&StdDevX);
  584. X        }
  585. X        if (AutoRange == FALSE) {
  586. X            StatisticsOfY(Data,NTuples,TupleSize,&MedianY,&StdDevY);
  587. X        }
  588. X        plaxes(
  589. X            MedianX,MedianY,
  590. X            Xopt, xtick, nxsub,
  591. X            Yopt, ytick,nysub
  592. X        );
  593. X    } else if (isInterval(V_Origin)) {
  594. X        x0 = V_Origin.val_u.udt_interval.int_lo;
  595. X        y0 = V_Origin.val_u.udt_interval.int_hi;
  596. X        plaxes(
  597. X            x0,y0,
  598. X            Xopt, xtick, nxsub,
  599. X            Yopt, ytick,nysub
  600. X        );
  601. X    } else if (isString(V_Origin) && Vstrcmp(V_Origin, "Automatic")) {
  602. X        plbox(Xopt, xtick, nxsub, Yopt, ytick,nysub);
  603. X    } else {
  604. X        plbox(Xopt, xtick, nxsub, Yopt, ytick,nysub);
  605. X    }
  606. X
  607. X    /* gridding    */
  608. X    if (VtoBoolean(V_Gridding)) {
  609. X        /* gridding on */
  610. X        plstyl(1, &mark1, &space1);
  611. X        plbox("g", xtick, nxsub, "g", ytick,nysub);
  612. X        plstyl(0, &mark0, &space0);
  613. X    }
  614. X
  615. X    /* plot curves    */
  616. X    if ((X = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
  617. X        perror("(Plot2D) ");
  618. X        ErrorExit();
  619. X    }
  620. X    if ((Y = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
  621. X        perror("(Plot2D) ");
  622. X        ErrorExit();
  623. X    }
  624. X    if (VtoBoolean(V_SupplyAbscissa)) {
  625. X        for (i=0; i<NTuples; i++)
  626. X            X[i] = LogX? log10((double)(i+1)) : (FLOAT)i;
  627. X    } else {
  628. X        for (i=0; i<NTuples; i++)
  629. X            X[i] = LogX? log10((double)Data[i][0]) : Data[i][0];
  630. X    }
  631. X    for (i=1; i<TupleSize; i++) {
  632. X        if (NLineStyles > 0) {
  633. X            /* use user-supplied patterns */
  634. X            plstyl(
  635. X                LineStyles[(i-1)%NLineStyles].ls_nelms,
  636. X                LineStyles[(i-1)%NLineStyles].ls_mark,
  637. X                LineStyles[(i-1)%NLineStyles].ls_space
  638. X            );
  639. X        } else {
  640. X            pllsty((i-1)%8 + 1);
  641. X        }
  642. X        plcol(i);
  643. X        for (j=0; j<NTuples; j++)
  644. X            Y[j] = LogY? log10((double)Data[j][i]) : Data[j][i];
  645. X
  646. X        if (VtoBoolean(V_PlotJoined))
  647. X            plline(NTuples, X, Y);
  648. X
  649. X        if (VtoBoolean(V_PlotPoints)) {
  650. X            plschr(0., VtoDbl(V_PointScale));
  651. X            pllsty(1);
  652. X            if (NPointCodes > 0) {
  653. X                plpoin(NTuples, X, Y, PointCodes[(i-1)%NPointCodes]);
  654. X            } else {
  655. X                plpoin(NTuples, X, Y, i);
  656. X            }
  657. X        }
  658. X    }
  659. X
  660. X
  661. X    /* restore line styles etc... */
  662. X    pllsty(1);
  663. X    plcol(1);
  664. X    plschr(0., VtoDbl(V_LabelScale));
  665. X    pllab(VtoString(V_XLabel), VtoString(V_YLabel), "");
  666. X    plschr(0., VtoDbl(V_TitleScale));
  667. X    pllab("","", VtoString(V_PlotTitle));
  668. X}
  669. X
  670. X#define    DTR(D)    (double)((D)*3.1415927/180.0)
  671. X
  672. Xvoid
  673. XPolarPlot(Data, NTuples, TupleSize)
  674. XFLOAT    **Data;
  675. Xint    NTuples;
  676. Xint    TupleSize;
  677. X{
  678. X    int    i, j;
  679. X    FLOAT    *A, *R;
  680. X    FLOAT    x, y;
  681. X    FLOAT    xtick, ytick;
  682. X    FLOAT    atick, rtick;
  683. X    FLOAT    Xmin,Xmax, Ymin,Ymax;
  684. X    FLOAT    r, rmax, ang, amax;
  685. X    FLOAT    MedianX, StdDevX;
  686. X    FLOAT    MedianY, StdDevY;
  687. X    FLOAT    da;
  688. X    FLOAT    x0,y0, xn,yn;
  689. X    double    pi;
  690. X    char    s[32];
  691. X    int    nxsub, nysub, nrsub, nasub;
  692. X    bool    InRadians, XisAngular;
  693. X    bool    AutoDomain, AutoRange;
  694. X#ifdef    ANSI_C
  695. X    void    StatisticsOfX(
  696. X            FLOAT **,
  697. X            int ,
  698. X            int ,
  699. X            FLOAT *,
  700. X            FLOAT *
  701. X        );
  702. X    void    StatisticsOfY(
  703. X            FLOAT **,
  704. X            int ,
  705. X            int ,
  706. X            FLOAT *,
  707. X            FLOAT *
  708. X        );
  709. X    double    log10(double); double fabs();
  710. X    double    sqrt(), sin(), cos(), atan();
  711. X#else
  712. X    void    StatisticsOfX();
  713. X    void    StatisticsOfY();
  714. X    double    log10(); double fabs();
  715. X    double    sqrt(), sin(), cos(), atan();
  716. X#endif
  717. X
  718. X    pi = 4.0 * atan((double)1.0);
  719. X    InRadians = Vstrcmp(V_AngularUnit, "degrees")? FALSE : TRUE;
  720. X    XisAngular = Vstrcmp(V_PolarVariable,"angle")? TRUE : FALSE;
  721. X
  722. X    pladv(0);
  723. X
  724. X    if (isString(V_AspectRatio)) {
  725. X        /* automatic --> STRETCH */
  726. X        plvasp(1.0);
  727. X    } else {
  728. X        assert(isDbl(V_AspectRatio));
  729. X        plvasp(VtoDbl(V_AspectRatio));
  730. X    }
  731. X
  732. X    /* window bounds */
  733. X    AutoDomain = FALSE;
  734. X    if (isString(V_Domain) && Vstrcmp(V_Domain, "All")) {
  735. X        /* all points    */
  736. X        Xmin = UserRect.XMin;
  737. X        Xmax = UserRect.XMax;
  738. X    } else if (isString(V_Domain) && Vstrcmp(V_Domain,"Automatic")) {
  739. X        /* automatic     */
  740. X        AutoDomain = TRUE;
  741. X        StatisticsOfX(Data, NTuples, TupleSize, &MedianX, &StdDevX);
  742. X        x = MedianX - 3.0*StdDevX;
  743. X        Xmin = UserRect.XMin < x? x : UserRect.XMin;
  744. X        x = MedianX + 3.0*StdDevX;
  745. X        Xmax = UserRect.XMax > x? x : UserRect.XMax;
  746. X    } else {
  747. X        assert(isInterval(V_Domain));
  748. X        Xmin = V_Domain.val_u.udt_interval.int_lo;
  749. X        Xmax = V_Domain.val_u.udt_interval.int_hi;
  750. X    }
  751. X
  752. X    AutoRange = FALSE;
  753. X    if (isString(V_Range) && Vstrcmp(V_Range, "All")) {
  754. X        /* all */
  755. X        Ymin = UserRect.YMin;
  756. X        Ymax = UserRect.YMax;
  757. X    } else if (isString(V_Range) && Vstrcmp(V_Range,"Automatic")) {
  758. X        /* automatic     */
  759. X        AutoRange = TRUE;
  760. X        StatisticsOfY(Data, NTuples, TupleSize, &MedianY, &StdDevY);
  761. X        y = MedianY - 3.0*StdDevY;
  762. X        Ymin = UserRect.YMin < y? y : UserRect.YMin;
  763. X        y = MedianY + 3.0*StdDevY;
  764. X        Ymax = UserRect.YMax > y? y : UserRect.YMax;
  765. X    } else {
  766. X        assert(isInterval(V_Range));
  767. X        Ymin = V_Range.val_u.udt_interval.int_lo;
  768. X        Ymax = V_Range.val_u.udt_interval.int_hi;
  769. X    }
  770. X    rmax = XisAngular? Ymax : Xmax;
  771. X    plwind(
  772. X        -1.3*rmax,
  773. X        1.3*rmax,
  774. X        -1.3*rmax,
  775. X        1.3*rmax
  776. X    );
  777. X
  778. X
  779. X    /* X and Y ticks */
  780. X    if (isString(V_XTick)) {
  781. X        /* automatic */
  782. X        xtick = XisAngular? (InRadians? pi/6.0 : 30.0) : 0.25 * rmax;
  783. X        nxsub = 4;
  784. X    } else {
  785. X        xtick = V_XTick.val_u.udt_interval.int_lo;
  786. X        nxsub = V_XTick.val_u.udt_interval.int_hi;
  787. X    }
  788. X    if (isString(V_YTick)) {
  789. X        /* automatic */
  790. X        ytick = XisAngular? 0.25 * rmax : (InRadians? pi/6.0 : 30.0);
  791. X        nysub = 4;
  792. X    } else {
  793. X        ytick = V_YTick.val_u.udt_interval.int_lo;
  794. X        nysub = V_YTick.val_u.udt_interval.int_hi;
  795. X    }
  796. X    rtick = XisAngular? ytick : xtick;
  797. X    atick = XisAngular? xtick : ytick;
  798. X    nasub = XisAngular? nxsub : nysub;
  799. X    nrsub = XisAngular? nysub : nxsub;
  800. X    amax = InRadians? 2.0 * pi : 360.0;
  801. X
  802. X    /* Axis options    */
  803. X    if (VtoBoolean(V_Gridding)) {
  804. X        /* draw circles    */
  805. X        plstyl(0, &mark1, space1);
  806. X        /*  circular ticks    */
  807. X        for (i=0,r=(rtick/nrsub); r<rmax; r += (rtick/nrsub),i++) {
  808. X            da = i%nrsub? 0.05 * atick : 0.1 * atick;
  809. X            for (ang=0; ang<amax; ang += atick) {
  810. X                x0 = r * cos(InRadians? ang-da : DTR(ang-da)); 
  811. X                y0 = r * sin(InRadians? ang-da : DTR(ang-da)); 
  812. X                xn = r * cos(InRadians? ang+da : DTR(ang+da)); 
  813. X                yn = r * sin(InRadians? ang+da : DTR(ang+da)); 
  814. X                pljoin(x0,y0,xn,yn);
  815. X            }
  816. X        }
  817. X#ifdef    CIRCULAR_TICKS
  818. X        /* major circles    */
  819. X        for (r=rtick; r<rmax; r += rtick) {
  820. X            x0 = r; y0 = 0.0;
  821. X            for (ang=1.0; ang<=360.0; ang++) {
  822. X                xn = r * cos((double)(DTR(ang))); 
  823. X                yn = r * sin((double)(DTR(ang))); 
  824. X                pljoin(x0,y0,xn,yn);
  825. X                x0 = xn; y0 = yn;
  826. X            }
  827. X        }
  828. X#endif
  829. X        /* draw spokes    */
  830. X        for (i=0,ang=0.0; ang<amax; ang += (atick/nasub),i++) {
  831. X            xn = rmax * cos(InRadians? ang : DTR(ang));
  832. X            yn = rmax * sin(InRadians? ang : DTR(ang));
  833. X            if (i%nasub == 0) {
  834. X                pljoin((FLOAT)0.0,(FLOAT)0.0,xn,yn);
  835. X            } else {
  836. X                pljoin(xn,yn,1.05*xn,1.05*yn);
  837. X            }
  838. X        }
  839. X    }
  840. X    /* write angle labels    */
  841. X    plschr(0., VtoDbl(V_AnnotationScale));
  842. X    j = amax / atick;
  843. X    for (ang=0.0,i=0; i<j; ang += atick,i++) {
  844. X        if (InRadians) {
  845. X            xn = rmax * cos((double)ang);
  846. X            yn = rmax * sin((double)ang);
  847. X            if (fabs(ang-(2.0*pi)) > 1.0e-2)
  848. X                sprintf(s, "%.2f #gp", (float)ang/pi);
  849. X        } else {
  850. X            xn = rmax * cos((double)DTR(ang));
  851. X            yn = rmax * sin((double)DTR(ang));
  852. X            sprintf(s, "%d", (int)(ang+0.5));
  853. X        }
  854. X        if (xn >= 0)
  855. X            plptex(xn,yn,xn,yn,-0.15, s);
  856. X        else
  857. X            plptex(xn,yn,-xn,-yn,1.15, s);
  858. X    }
  859. X    plstyl(0, &mark0, &space0);
  860. X
  861. X    /* plot curves    */
  862. X    if ((A = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
  863. X        perror("(Plot2D) ");
  864. X        ErrorExit();
  865. X    }
  866. X    if ((R = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
  867. X        perror("(Plot2D) ");
  868. X        ErrorExit();
  869. X    }
  870. X    if (VtoBoolean(V_SupplyAbscissa)) {
  871. X        for (i=0; i<NTuples; i++)
  872. X            if (XisAngular)
  873. X                A[i] = i*2.0*pi/(NTuples-1);
  874. X            else
  875. X                R[i] = i*rmax/(NTuples-1);
  876. X    } else {
  877. X        for (i=0; i<NTuples; i++)
  878. X            if (XisAngular)
  879. X                A[i] = InRadians? Data[i][0] :  DTR(Data[i][0]);
  880. X            else
  881. X                R[i] = Data[i][0];
  882. X    }
  883. X    for (i=1; i<TupleSize; i++) {
  884. X        if (NLineStyles > 0) {
  885. X            /* use user-supplied patterns */
  886. X            plstyl(
  887. X                LineStyles[(i-1)%NLineStyles].ls_nelms,
  888. X                LineStyles[(i-1)%NLineStyles].ls_mark,
  889. X                LineStyles[(i-1)%NLineStyles].ls_space
  890. X            );
  891. X        } else {
  892. X            pllsty((i-1)%8);
  893. X        }
  894. X        plcol(i);
  895. X
  896. X        /* get and convert data if necessary    */
  897. X        if (InRadians && !XisAngular)
  898. X            for (j=0; j<NTuples; j++)
  899. X                A[j] = Data[j][i];
  900. X        else if ((!InRadians) && !XisAngular)
  901. X            for (j=0; j<NTuples; j++)
  902. X                A[j] = DTR(Data[j][i]);
  903. X        else
  904. X            for (j=0; j<NTuples; j++)
  905. X                R[j] = Data[j][i];
  906. X
  907. X        /* convert to cartesian    A<->X  R<->Y*/
  908. X        for (j=0; j<NTuples; j++) {
  909. X            x0 = R[j] * cos(A[j]);
  910. X            y0 = R[j] * sin(A[j]);
  911. X            A[j] = x0;
  912. X            R[j] = y0;
  913. X        }
  914. X        if (VtoBoolean(V_PlotJoined)) {
  915. X            x0 = A[0];
  916. X            y0 = R[0];
  917. X            for (j=1; j<NTuples; j++) {
  918. X                pljoin(A[j-1], R[j-1], A[j],R[j]);
  919. X            }
  920. X        }
  921. X
  922. X        if (VtoBoolean(V_PlotPoints)) {
  923. X            plschr(0., VtoDbl(V_PointScale));
  924. X            pllsty(1);
  925. X            if (NPointCodes > 0) {
  926. X                plpoin(NTuples, A, R, PointCodes[(i-1)%NPointCodes]);
  927. X            } else {
  928. X                plpoin(NTuples, A, R, i);
  929. X            }
  930. X        }
  931. X    }
  932. X
  933. X    /* restore line styles etc... */
  934. X    pllsty(1);
  935. X    plcol(1);
  936. X    plschr(0., VtoDbl(V_LabelScale));
  937. X    pllab(VtoString(V_XLabel), VtoString(V_YLabel), "");
  938. X    plschr(0., VtoDbl(V_TitleScale));
  939. X    pllab("","", VtoString(V_PlotTitle));
  940. X}
  941. X
  942. X
  943. Xvoid
  944. XStatisticsOfX(Data, NTuples,TupleSize, Median, StdDev)
  945. XFLOAT    **Data;
  946. Xint    NTuples;
  947. Xint    TupleSize;
  948. XFLOAT    *Median;
  949. XFLOAT    *StdDev;
  950. X{
  951. X#ifdef    ANSI_C
  952. X    FLOAT    MedianOfX();
  953. X    FLOAT    StdDevOfX();
  954. X/*debug
  955. X    FLOAT    MedianOfX(FLOAT **Data,int NTuples,int TupleSize);
  956. X    FLOAT    StdDevOfX(FLOAT **Data,int NTuples,int TupleSize);
  957. Xdebug*/
  958. X#else
  959. X    FLOAT    MedianOfX();
  960. X    FLOAT    StdDevOfX();
  961. X#endif
  962. X    *StdDev = StdDevOfX(Data,NTuples,TupleSize);    
  963. X    *Median = MedianOfX(Data,NTuples,TupleSize);
  964. X}
  965. X
  966. X
  967. Xvoid
  968. XStatisticsOfY(Data, NTuples,TupleSize, Median, StdDev)
  969. XFLOAT    **Data;
  970. Xint    NTuples;
  971. Xint    TupleSize;
  972. XFLOAT    *Median;
  973. XFLOAT    *StdDev;
  974. X{
  975. X#ifdef    ANSI_C
  976. X    FLOAT    MedianOfY();
  977. X    FLOAT    StdDevOfY();
  978. X/*debug
  979. X    FLOAT    MedianOfY(FLOAT **Data,int NTuples,int TupleSize);
  980. X    FLOAT    StdDevOfY(FLOAT **Data,int NTuples,int TupleSize);
  981. Xdebug*/
  982. X#else
  983. X    FLOAT    MedianOfY();
  984. X    FLOAT    StdDevOfY();
  985. X#endif
  986. X    *StdDev = StdDevOfY(Data,NTuples,TupleSize);    
  987. X    *Median = MedianOfY(Data,NTuples,TupleSize);
  988. X}
  989. X
  990. X
  991. XFLOAT
  992. XStdDevOfX(Data,NTuples,TupleSize)
  993. XFLOAT    **Data;
  994. Xint    NTuples;
  995. Xint    TupleSize;
  996. X{
  997. X    int    i;
  998. X    FLOAT    StdDev, x;
  999. X    FLOAT    var, sum, mean;
  1000. X    double sqrt(), fabs();
  1001. X
  1002. X    if (VtoBoolean(V_SupplyAbscissa)) {
  1003. X        return((FLOAT)(NTuples*NTuples)/12.0);
  1004. X    }
  1005. X
  1006. X    sum = 0.0;
  1007. X    for (i=0; i<NTuples; i++)
  1008. X        sum += Data[i][0];
  1009. X    mean = sum / NTuples;
  1010. X
  1011. X    var = 0.0;
  1012. X    for (i=0; i<NTuples; i++) {
  1013. X        x = fabs((double)(Data[i][0] -  mean));
  1014. X        var += x * x;
  1015. X    }
  1016. X    var /= (NTuples-1);
  1017. X    StdDev = sqrt((double)var);
  1018. X    return(StdDev);    
  1019. X}
  1020. X
  1021. X
  1022. XFLOAT
  1023. XStdDevOfY(Data,NTuples,TupleSize)
  1024. XFLOAT    **Data;
  1025. Xint    NTuples;
  1026. Xint    TupleSize;
  1027. X{
  1028. X    int    i,j,D;
  1029. X    FLOAT    StdDev, x;
  1030. X    FLOAT    var, sum, mean;
  1031. X    FLOAT    n;
  1032. X    double sqrt(), fabs();
  1033. X
  1034. X    D = VtoBoolean(V_SupplyAbscissa) ? 0 : 1;
  1035. X    n = NTuples * (TupleSize-D);
  1036. X    sum = 0.0;
  1037. X    for (i=0; i<NTuples; i++)
  1038. X        for (j=D; j<TupleSize; j++) {
  1039. X            sum += Data[i][j];
  1040. X        }
  1041. X    mean = sum / n;
  1042. X
  1043. X    var = 0.0;
  1044. X    for (i=0; i<NTuples; i++)
  1045. X        for (j=D; j<TupleSize; j++) {
  1046. X            x = fabs((double)(Data[i][j] -  mean));
  1047. X            var += x * x;
  1048. X        }
  1049. X    var /= (n-1);
  1050. X    StdDev = sqrt((double)var);
  1051. X    return(StdDev);    
  1052. X}
  1053. X
  1054. X#define    BIG    1.0e30
  1055. X#define    AFAC    1.5
  1056. X#define    AMP    1.5
  1057. X
  1058. XFLOAT
  1059. XMedianOfX(Data,NTuples,TupleSize)
  1060. XFLOAT    **Data;
  1061. Xint    NTuples;
  1062. Xint    TupleSize;
  1063. X/* an iterative computation of the median...
  1064. X * Code taken from Numerical Recipes..
  1065. X */
  1066. X{
  1067. X    int    np, nm, i;
  1068. X    FLOAT    xx,xp,xm,sumx,sum,eps,stemp,dum,ap,am,aa,a;
  1069. X    FLOAT    Median;
  1070. X    double    fabs();
  1071. X
  1072. X    if (VtoBoolean(V_SupplyAbscissa))
  1073. X        return((FLOAT)(0.5*NTuples));
  1074. X
  1075. X    /* first estimate    */
  1076. X    a = 0.5 * (Data[0][0] + Data[NTuples-1][0]);
  1077. X    eps = fabs((double)(Data[NTuples-1][0] - Data[0][0]));
  1078. X
  1079. X    am = -(ap=BIG);
  1080. X    for (;;) {
  1081. X        sum = sumx = 0.0;
  1082. X        np = nm = 0;        /* # point above and below median */
  1083. X        xm = -(xp = BIG);
  1084. X
  1085. X        for (i=0; i<NTuples; i++) {
  1086. X            xx = Data[i][0];
  1087. X            if (xx != a) {
  1088. X                if (xx > a) {
  1089. X                    ++np;
  1090. X                    if (xx < xp) xp = xx;
  1091. X                } else if (xx < a) {
  1092. X                    ++nm;
  1093. X                    if (xx > xm) xm = xx;
  1094. X                }
  1095. X                sum += dum=1.0/(eps+fabs((double)(xx-a)));
  1096. X                sumx += xx * dum;
  1097. X            }
  1098. X        }
  1099. X
  1100. X        stemp = (sumx / sum) - a;
  1101. X        if ((np-nm) >= 2) {
  1102. X            /* guess is too low make another pass */
  1103. X            am = a;
  1104. X            aa = stemp < 0.0? xp : xp + stemp*AMP;
  1105. X            if (aa > ap) aa = 0.5*(a+ap);
  1106. X            eps = AFAC * fabs((double)(aa-a));
  1107. X            a = aa;
  1108. X        } else if ((nm-np) >= 2) {
  1109. X            /* guess is too high make another pass */
  1110. X            ap = a;
  1111. X            aa = stemp > 0.0? xm : xm+stemp*AMP;
  1112. X            if (aa < am) aa = 0.5*(a+am);
  1113. X            eps = AFAC*fabs((double)(aa-a));
  1114. X            a = aa;
  1115. X        } else {    /* got it    */
  1116. X            if (NTuples%2 == 0) {
  1117. X                Median = 0.5*(np==nm? xp+xm : np>nm? a+xp : xm+a);
  1118. X            } else {
  1119. X                Median = np == nm? a : np>nm? xp : xm;
  1120. X            }
  1121. X            return(Median);
  1122. X        }
  1123. X    }
  1124. X}
  1125. X
  1126. XFLOAT
  1127. XMedianOfY(Data,N,M)
  1128. XFLOAT    **Data;
  1129. Xint    N;
  1130. Xint    M;
  1131. X/* an iterative computation of the median...
  1132. X * Code taken from Numerical Recipes..
  1133. X */
  1134. X{
  1135. X    int    n, np, nm, i,j,D;
  1136. X    FLOAT    xx,xp,xm,sumx,sum,eps,stemp,dum,ap,am,aa,a;
  1137. X    FLOAT    Median;
  1138. X    double    fabs();
  1139. X
  1140. X    D = VtoBoolean(V_SupplyAbscissa) ? 0 : 1;
  1141. X    n = N * (M-D);
  1142. X
  1143. X    /* first estimate    */
  1144. X    a = 0.5 * (Data[0][D] + Data[N-1][M-1]);
  1145. X    eps = fabs((double)(Data[N-1][M-1] - Data[0][D]));
  1146. X
  1147. X    am = -(ap=BIG);
  1148. X    for (;;) {
  1149. X        sum = sumx = 0.0;
  1150. X        np = nm = 0;        /* # point above and below median */
  1151. X        xm = -(xp = BIG);
  1152. X
  1153. X        for (i=0; i<N; i++) {
  1154. X            for (j=D; j<M; j++) {
  1155. X                xx = Data[i][j];
  1156. X                if (xx != a) {
  1157. X                    if (xx > a) {
  1158. X                        ++np;
  1159. X                        if (xx < xp) xp = xx;
  1160. X                    } else if (xx < a) {
  1161. X                        ++nm;
  1162. X                        if (xx > xm) xm = xx;
  1163. X                    }
  1164. X                    sum += dum=1.0/(eps+fabs((double)(xx-a)));
  1165. X                    sumx += xx * dum;
  1166. X                }
  1167. X            }
  1168. X        }
  1169. X
  1170. X        stemp = (sumx / sum) - a;
  1171. X        if ((np-nm) >= 2) {
  1172. X            /* guess is too low make another pass */
  1173. X            am = a;
  1174. X            aa = stemp < 0.0? xp : xp + stemp*AMP;
  1175. X            if (aa > ap) aa = 0.5*(a+ap);
  1176. X            eps = AFAC * fabs((double)(aa-a));
  1177. X            a = aa;
  1178. X        } else if ((nm-np) >= 2) {
  1179. X            /* guess is too high make another pass */
  1180. X            ap = a;
  1181. X            aa = stemp > 0.0? xm : xm+stemp*AMP;
  1182. X            if (aa < am) aa = 0.5*(a+am);
  1183. X            eps = AFAC*fabs((double)(aa-a));
  1184. X            a = aa;
  1185. X        } else {    /* got it    */
  1186. X            if (n%2 == 0) {
  1187. X                Median = 0.5*(np==nm? xp+xm : np>nm? a+xp : xm+a);
  1188. X            } else {
  1189. X                Median = np == nm? a : np>nm? xp : xm;
  1190. X            }
  1191. X            return(Median);
  1192. X        }
  1193. X    }
  1194. X}
  1195. END_OF_FILE
  1196. if test 25452 -ne `wc -c <'Csrc/plotting.c'`; then
  1197.     echo shar: \"'Csrc/plotting.c'\" unpacked with wrong size!
  1198. fi
  1199. chmod +x 'Csrc/plotting.c'
  1200. # end of 'Csrc/plotting.c'
  1201. fi
  1202. echo shar: End of archive 3 \(of 3\).
  1203. cp /dev/null ark3isdone
  1204. MISSING=""
  1205. for I in 1 2 3 ; do
  1206.     if test ! -f ark${I}isdone ; then
  1207.     MISSING="${MISSING} ${I}"
  1208.     fi
  1209. done
  1210. if test "${MISSING}" = "" ; then
  1211.     echo You have unpacked all 3 archives.
  1212.     rm -f ark[1-9]isdone
  1213. else
  1214.     echo You still need to unpack the following archives:
  1215.     echo "        " ${MISSING}
  1216. fi
  1217. ##  End of shell archive.
  1218. exit 0
  1219. -- 
  1220. Mail submissions (sources or binaries) to <amiga@uunet.uu.net>.
  1221. Mail comments to the moderator at <amiga-request@uunet.uu.net>.
  1222. Post requests for sources, and general discussion to comp.sys.amiga.
  1223.